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Abstract 

The 1-D, quasi 1-D and 2-D Euler solvers based 
on the method of space-time conservation element 
and solution element are used to simulate various 
flow phenomena including shock waves, Mach stem, 
contact surface, expansion waves, and their inter- 
sections and reflections. Seven test problems are 
solved to demonstrate the capability of this method 
for handling unsteady compressible flows in various 
configurations. Numerical results so obtained are 
compared with exact solutions and/or numerical so- 
lutions obtained by schemes based on other estab- 
lished computational techniques. Comparisons show 
that the present Euler solvers can generate highly 
accurate numerical solutions to complex flow prob- 
lems in a straightforw^ard manner without using any 
ad hoc techniques in the scheme. 

1. Introduction 

The method of space-time conservation element 
and solution element (to be abbreviated as the 
CE/SE method) is a new numerical method devel- 
oped by Chang for solving conservation laws.^”'^ It 
is different in both concept and methodology from 
the well-established traditional methods such as the 
finite difference, finite volume, finite element and 
spectral methods. It is designed from a physicist’s 
perspective to overcome several key limitations of 
the traditional methods. 

Simplicity, generality and accuracy are weighted in 
the development of the present method while satis- 
fying the fundamental computational requirements. 
Its salient properties are summarized briefly as fol- 
lows. First, the concepts of conservation element 
and solution element are introduced to enforce both 
local and global flux conservation in space and time 
instead of in space only. Second, all the dependent 
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variables and their spatial derivatives are consid- 
ered as individual unknowns to be solved simultane- 
ously at each grid point. Third, no approximation 
techniques other than Taylor’s series expansion, no 
monotonous constraints, and no characteristic-based 
techniques are used in the design of the scheme. De- 
tailed description of this method and the accompa- 
nied analysis are referred to Refs. 1 and 2. 

Various efficient numerical schemes based on the 
CE/SE method have been developed for solving dif- 
ferent flow problems, especially the problems in the 
presence of shock waves with discontinuous flow 
properties. Of those schemes, the one- and two- 
dimensional time marching Euler solvers are em- 
ployed here to solve problems involving flow phe- 
nomena that are more complex than those shown in 
Refs. 1 and 2. In addition, a quasi one-dimensional 
Euler solver is constructed in this work which is 
aimed at dealing with problems in axisymmetric con- 
figurations. Comparisons of the computed results 
with published data are made to demonstrate the 
simplicity and accuracy of the present Euler solvers. 

Detailed description of the governing equations 
and CE/SE 1-D, quasi 1-D, 2-D Euler solvers used 
in the following numerical tests is referred to Ref. 3, 
while numerical results and discussions are described 
in details as follows. 

2. Numerical Results 
2.1) 1-D shock-tube problems 

In the 1-D weighted-average Euler a-e scheme used 
in the numerical tests for the following four shock- 
tube problems, the parameter e is set as 0.5, while 
a is 1 except in the first problem where tv is 4. 

a) The Lax problem 

The initial conditions in the region [—8,6] on the 
X axis are defined as 

Pi = 0.445, ui = 0.6989, pi = 3.5277 x < 0 (1) 
p,. = 0.5, = 0.0, Pr = 0.571 X > 0 (2) 
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The numerical solution at i = lOOAt obtained by 
using the present scheme, under the same compu- 
tational conditions (CFL=0.95 and Ax = 0.1) as 
those used in Ref. 5, is shown with the exact solu- 
tion represented by solid lines in Fig. 1. Compar- 
isons were made by Harten^ to appraise four numeri- 
cal schemes, namely, the ROE, LW (Lax-Wendroff), 
ULTl, and ULTIC schemes, which were used to 
solve the same problem. The last two are the TVNI 
finite-difference schemes designed by Harten. The 
numerical results plotted in Figs. 2(a)-(d) of Ref. 5 
are not reproduced here. Among the four solutions, 
the LW solution is the worst showing serious oscil- 
lations near the shock discontinuity. In the ROE 
solution, the oscillations disappear but up to four 
grid points are needed to resolve a shock wave, while 
a significant smear is found either near the contact 
surface or within the expansion fan. Much improve- 
ment is revealed in the ULTl solution, which re- 
solves a shock discontinuity with only 2 grid points 
and shows an excellent agreement with the exact so- 
lution except in the diffusive region near the contact 
surface. The ULTl solution closely resembles the 
solution plotted in Fig, 1 obtained by use of the 
present Euler solver. Harten showed that a further 
refinement can be obtained by using the higher order 
accurate ULTIC scheme. 

b) The Sjogreen problem 

This problem is taken from Ref. 6, whose initial 
conditions are 

/?/ = 1.0, = p, = 0.4 0<x<0.5 (3) 

Pr - TO, Ur = 2, Pr - 0.4 0.5 < X < 1.0 (4) 

The initial velocity discontinuity causes two rar- 
efaction waves to propagate in opposite directions, 
leaving in between a region of high vacuum. It 
was mentioned^ that several Godunov-type schemes 
failed in this problem due to the extremely low pres- 
sure in the middle region. The CE/SE solution at 
t = bOAt based on 100 grid points and At = 0.002 
is shown in Fig. 2, in which the exact solution is 
represented by solid lines. It can be seen that the 
present solution agrees very well with the exact so- 
lution, without showing negative pressure values in 
the middle region. The solution displays an accu- 
racy which is comparable to that obtained by Xu et 
al.^ using a gas-kinetic scheme with 200 grid points. 

c) The Shu-Osher problem 

Examined in this problem is the interaction of a 
moving shock of Ms — 3 with a sinusoidal density 
wave.® The initial conditions in the region [—5, 5] are 
described as 

= 3.857, = 2.629, p/ = 10.333 x < ~4 (5) 


= 1 -f 0.2 sin 5x, Ur = 0, Pr ~ I otherwise (6) 

This problem does not have an exact solution. 
Several upwind schemes have been used to solve this 
problem to compare their abilities in resolving the 
peaks appeared in the solution.^ The CE/SE solu- 
tion at t=1.8 obtained by using 800 grid points with 
At — 0.0015 (CFL=0.582) is shown in Fig. 3. The 
present solution is comparable to those obtained in 
Ref. 9 by using the TVDl and TVD2 schemes with 
the same number of grid points. 

d) The Woodward-Colella problem 

This problem, concerning the interaction of two 
blast waves in a close-ended tube, was proposed by 
Woodward and Colella without an exact solution. 
The initial conditions are 

Pi — 1.0, uj =0, Pi = 1000 X < 0.1 (7) 

Pm — 1-0, Um — 0, Pm = 0.01 0.1 < X < 0.9 (8) 

pr ~ 1.0, Uy =0, Pr — 100 0-9 < J' < 1.0 (9) 

The two ends are at x — 0 and x \ where t he re- 
flecting boundary conditions are imposed. Detailed 
reflecting boundary conditions used in the present 
scheme are referred to Ref. 4. The CE/SE so- 
lution at t=0.038 based on 800 grid points with 
At = 1.25 X 10-^(CFL=0.3524) is shown in Fig. 4. 
The flowfield at t=0.038 contains three contact sur- 
faces and two shock waves. It can be seen that the 
contact surfaces are much smeared than the shock 
discontinuities. A comparison with the numerical 
solutions obtained by using AUSM+, Roe, Van leer, 
AUSMDV and AUSM^-w splitting schemes for the 
same problem shown in Fig, 5 of Ref. 11 reveals that 
the CE/SE solution is at least of the same accuracy. 

It hcLs been demonstrated that the present 1-D Eu- 
ler solver can generate highly accurate solutions to 
shock-tube problems involving various discontinu- 
ities, even though it does not need the implement 
of monotonous restraints, TVNI, and entropy con- 
ditions cLs did in Ref. 5. This simple scheme can 
be used without difficulty to solve any 1-D problems 
governed by Euler equations. 

2.2) A quasi 1-D nozzle flow 

An axisymmctric nozzle with cross-sectional area 
A(x) = 1.398 -h 0.347 tanh(0. 8a; — 4) described in 
Refs. 5 and 12 is reconsidered here. Numerical solu- 
tions obtained by use of the qucisi 1-D Euler solver 
for CFL=0.9 with 20 and 32 grid points are shown 
in Figs. 5 and 6, respectively. The present solution 
with 20 grid points is at least as good as those ob- 
tained by ROE and ULTl schemes shown in Fig. 6 
of Ref. 5 with tlie same number of grid points, while 
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the present solution with 32 grid points is better 
than the ULTl solution with 50 grid points. 

2.3) 2-D supersonic flow past a step 

Consider the supersonic channel flow past a step 
depicted in Fig. 7(a). The flow exhibits complicated 
phenomena which include Mach stem, slip surface, 
shock wave, expansion fan, and their interactions 
and reflections. This is a standard benchmark prob- 
lem in the literature. It was used to test Marten’s 
TVNI ULTIC scheme,^ Giannakouros and Karni- 
adakis’ spectral element-FCT method, and Van 
Leer’s ultimate conservative difference scheme. It 
was also adopted by Woodward and Colella^^ to 
compare the accuracy of different numerical meth- 
ods in handling a shock discontinuity. This problem 
is solved again here to demonstrate the robustness 
of the CE/SE method. 

The 2-D weighted-average Euler a-e scheme with 
c = 0.5 and a = 1 is used in this numerical test. The 
grid distribution shown in Fig. 7(b) indicates that 
no grid point is placed at the upper corner of t.lie 
step to avoid the singular flow behavior there. The 
freestream condition is set at the inlet, while the 
condition at the exit is extrapolated from the inte- 
rior, The reflective boundary condition is imposed 
on solid walls. Detailed description of boundary con- 
ditions can be found in Ref. 3. 

The density contours in the solutions obtained by 
the present Euler solver with 61x21, 121x41 and 
241x81 grid points are shown in Fig. 8. Similar 
contour plots are displayed in Figs. 7(a)-(g) of Ref, 
10 based on six selected numerical schemes. Accord- 
ing to Ref. 10, the ranking of these six methods in 
terms of accuracy is as follows: PPM(both PPMLR 
and PPMDE), MUSCL, ETBFCT, BBC, MacCor- 
mack’s scheme, and Godunov’s scheme. 

Comparisons under the same computational con- 
ditions (CFL=0,8, At = 0.0025 for 241x81 grid) 
show that the present solution is as good as those 
obtained by the accurate MUSCL and PPMDE 
schemes and is much better than Godunov’s solu- 
tion. A direct comparison cannot be made with 
other mentioned schemes because of their different 
time step sizes and CFL numbers. The Mach stem, 
triple point, slip surface, expansion fan at the cor- 
ner, and the interaction between the reflected shock 
with rarefaction waves are accurately simulated in 
the present numerical solutions. 

2.4) A 2-D blast flowfield 

Considered here is a blast flowfield generated by 
an open-ended cylindrical shock tube, which was 
simulated in Ref. 15 using a TVD finite volume 


method with numerical techniques for controlling ar- 
tificial compressibility and dissipation. The flowfield 
involves complicated phenomena including vortex, 
blast wave, rarefaction wave, normal shock and their 
mutual interactions. The early time development of 
vortex and shock diffraction and the subsequent flow 
evolution were simulated in Ref. 15 up to 1.5 msec 
to show a fair comparison with experimental data. 
This problem is solved here on cartesian coordinates 
using the CE/SE method to demonstrate its versa- 
tility. 

The two-dimensional shock tube configuration 
adopted for numerical computation is depicted in 
Fig. 9, in which the blank space is used to represent 
the solid tube wall above the plane of symmetry on 
the X axis. Displayed schematically in the figure are 
some representative grid points. A shock wave is 
created by the sudden removal of a diaphragm at 
the lip of the tube which separates a compressed 
fluid in region 2 inside the tube from the surround- 
ing stagnant fluid in region 1. The initial conditions 
are described by 

Pi = 1/1.4, Pi = 1.0, ui — 0.0, v\ — 0.0 (10) 

P2 = 2.443, p2 = 2.28, «2 = 0.982, U2 = 0.0 (11) 

The 2-D Euler solver with e = 0.5 and a = 1 
is used to compute the flowfield with At = 0.0025 
on a mesh of 49x97 grid points. The nonreflective 
boundary condition is set at the inlet, outlet, and 
upper boundary of the computational domain, while 
the reflective boundary condition is set on the x axis 
and at all tube walls. 

Qualitative agreements between computed and 
measured positions of the vortex center (estimated 
from Fig. 15 of Ref. 15) at six time levels are 
revealed in Fig. 10, A quantitative comparison is 
not feasible due to the fact that the experiment was 
performed in an axisymmetric configuration whereas 
the computation was in 2-D. However, at an early 
stage the measured vortex positions are correctly 
predicted by the 2-D code. 

The velocity vector fields at six time levels are 
plotted in Fig. 11 to show the formation and move- 
ment of both the vortex and blast wave. The vortex 
can be recognized by its characteristic revolving flow 
pattern and the blast wave is represented by a shal- 
low banded curve. The computed flow features can 
be detected in photographs taken at different time 
steps. Our numerical results with 4753 grid points 
uniformly distributed in the computational domain 
are comparable to those reported in Ref. 15 bcised on 
a finite volume method containing 7377 cells, with 
much smaller cells densely distributed in the neigh- 
borhood of the tube exit. 
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After this test, the 2-D blast flowfield generated 
by the same initial conditions is simulated to further 
test the robustness of the present Euler solver. The 
computational domain is enlarged to — 1 < x < 3 
and 0 < y < 3, which is the same as that of the 
axisymmetric flowfield simulated in Ref. 15. The 
boundary conditions are the same as those described 
previously except the upper surface is now replaced 
by a solid wall. The reflection of the blast wave from 
the upper surface causes additional flow phenomena 
that were not observed in the experiment. 

To show the time history of flow development, nu- 
merical solutions at eight time levels obtained us- 
ing 161x121 grid points are shown in Figs. 12 and 
13, in which the pressure and density contours rang- 
ing from 0 to 5.88 with a constant interval of 0.049 
are plotted. The sequential plots reveal that as the 
blast wave initiated from the open end of the shock 
tube propagates away, a vortex is developed at the 
lip of the tube wall, which moves downstream with 
an ascending motion. When the blast wave reaches 
the upper wall, it is reflected as shown in the plot 
at t^l.O msec. In the meantime on the tube axis a 
normal shock is formed ahead the vortex and is mov- 
ing slowly in the downstream direction, while a jet 
shear layer is created at the lip of the shock tube. At 
t=:L5 msec, the portion of the blast wave that is re- 
flected from the upper wall is shown to move toward 
the vortex. After passing the vortex, the blast wave 
becomes curved and keeps moving forward to inter- 
act with the normal shock as shown in the plots at 
t=1.7 msec. At t=1.9 msec, the flow pattern reveals 
that as a result of the interaction, the blast wave is 
broken into two parts while several new vortices are 
created. More complex flow patterns are shown at 
t=2.1 and 2.3 msec, describing further reflection and 
interaction of shock waves arid vortices. 

Despite the difference between 2-D and axisym- 
metric configurations, the computed flow fields agree 
extremely well with those shown in the shadowgraph 
pictures of Ref. 15 at early time steps of t=0.1996 
and 0.4937 msec. An axisymmetric Euler solver 
based on the CE/SE method is being developed, 
which will be used to simulate the realistic exper- 
imental flow conditions. 

3. Conclusions 

The 1-D, quasi TD and 2-D Euler solvers have 
been validated using test problems. Numerical ex- 
amples have been used to compare the CE/SE Euler 
solvers with established schemes such as the TVNI 
schemes designed by Harten, the upwind schemes 
used by Woodward and Colella, and others. It has 
been demonstrated that the present. Euler solvers 
can generate highly accurate numerical solutions 


without requiring any special treatments for flow dis- 
continuities, such as the inclusion of artificial viscos- 
ity, blending of low- and high-order-accurate fluxes, 
the use of nonlinear solution to Riemann^s prob- 
lem as suggested in Ref. 10, or the TVD property 
used in Ref. 15. Its inherent features of simplicity, 
generality and accuracy indicate that, with further 
improvements, the space- time conservation element 
and solution element method may be developed to 
become a versatile tool for solving general fluid dy- 
namic problems. 
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Fig. 1 The CE/SE toluilon for Lax problem 
(crt-o.os, Ax-at. a-4. i»iooAt) 
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Fig. 2 The CE/SE eolution for SJogreen problem 
(At-O.OOZ. AX-0.01. e«l. l-60At) 
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Fig. 3 The CE/SE solution for Shu^Osher problem 
(M-O.OOtfl. AX-0.0120, a-1. I- 1.6) 
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rif. 4 The CE/SE eoIuUon for bUet weve problem "«■ * eteedp-eUU eoluUon of noxxle probelm 
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Fig. 8 Density contours for flow over a step at t=4 
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Fig. 12 Pressure contours at eight time levels 









Fig. 13 Density contours at eight time levels 









